Introducción:
La diferenciación celular es un proceso en el que intervienen varios factores a nivel genético y epigenético, que provocan diferencias en la regulación de la transcripción de grupos de genes distintos. Estas diferencias son cruciales para desarrollar funciones específicas en cada célula.
Entre los factores que influyen se encuentran la estructura de la cromatina y las modificaciones de histonas. Por ejemplo, la modificación H3K4me1 está presente en regiones con enhancers activos que ayudan a activar la transcripción de genes, mientras que la histona H3K4me3 es una marca relacionada con promotores activos.
Usando los datos de un experimeto de CHiP-Seq de la histona H3K4me1 obtenidos del proyecto BLUEPRINT, se busca obtener los promotores y enhancers asociados a esta histona en la célula, usando los datos obtenidos en el artículo “Lineage-Specific Genome Architecture Links Enhancers and Non-coding Disease Variants to Target Gene Promoters” de Javierre et al. (2016), en donde se utilizó el método de promotor capture Hi-C (PCHi-C) para estudiar la interacción espacial de la cromatina con promotores de genes específicos, a los que se refieren como “baits” y los elementos reguladores distales, como los enhancers, que interactúan con estos promotores, denominados “other ends”.
De acuerdo al artículo se espera obtener mayor enriquecimiento de enhacers en neutrófilos maduros. A partir de métodos estadísticos se pretende comprobar que, efectivamente, existe una diferencia sinificativa entre las variables categóricas “enhancer” y “promotores” con dicha histona asociada.
Finalmente, se analiza el enriquecimiento de los factores de transcripción que se unen en estas regiones. Esto a partir de matrices de peso obtenidas de las bases de JASPAR y usando secuencias permutadas como control.
Razón biológica de la marca H3k4me1
La H3K4me1 es una marca epigenética de metilación de la histona H3 en la posición 4 del residuo de lisina. Algunas características de la marca H3K4me1 son:
Marcador de regiones de regulación génica: La H3K4me1 es una marca epigenética comúnmente encontrada en regiones no codificantes del genoma, como los enhancers, que son regiones del ADN que regulan la expresión génica.
Asociada con la activación de la transcripción: La presencia de la marca H3K4me1 en enhancers y otras regiones reguladoras ha sido asociada con la activación de la transcripción génica, lo que sugiere que puede tener un papel importante en la regulación de la expresión génica.
Potencial biomarcador: La presencia de la marca H3K4me1 ha sido propuesta como un biomarcador para varios tipos de cáncer, lo que sugiere que puede tener un papel importante en la progresión del cáncer y la respuesta a la terapia.
En articulo, en la figura 3B, se observa el enriquecimiento de Promoting Interacting Regions (PIRs) asociadas con marcas de histona. En esta encontramos que la histona H3K4me1 se encuentra enriquecida en todos los tipos celulares.
Aunado a esto, en la figura 3D los neutrófilos muestran un mayor enriquecimiento de promotores y enhancers asociados.
Dado el potencial regulatorio de la histona, elegimos estudiarla en un neutrófilo para esperar observar un mayor enriquecimiento de enhancers.
Descripción de los datos
Interacciones promotor-enhancer
Tenemos una tabla delimitada por tabulaciones que enumera todas las interacciones (puntuación CHiCAGO >=5) detectadas entre promotores activos y enhancers (definidos sobre la base de Regulatory Build y las anotaciones de cromatina BLUEPRINT). El puntaje CHiCAGO es una medida estadística utilizada para evaluar la confiabilidad de los datos de interacción de cromatina obtenidos a partir de experimentos de HI-C.
interacciones <- read.table("ActivePromoterEnhancerLinks.tsv", header = T, sep = "\t")
head(interacciones)## baitChr baitSt baitEnd baitID oeChr oeSt oeEnd oeID
## 1 chr1 1206873 1212438 254 chr1 943676 957199 228
## 2 chr1 1206873 1212438 254 chr1 1034268 1040208 235
## 3 chr1 1206873 1212438 254 chr1 1040208 1043143 236
## 4 chr1 1206873 1212438 254 chr1 1069045 1083958 242
## 5 chr1 1206873 1212438 254 chr1 1083958 1091234 243
## 6 chr1 1206873 1212438 254 chr1 1585571 1619752 304
## cellType.s.
## 1 nCD8
## 2 nCD4,nCD8,Mac0,Mac1,Mac2,MK,Mon
## 3 nCD4,nCD8,Mac0,Mac1,Mac2,MK
## 4 nCD8
## 5 nCD8
## 6 Neu
## sample.s.
## 1 C0066PH1
## 2 S007DDH2,S007G7H4,C0066PH1,S00C2FH1,S00390H1,S001MJH1,S001S7H2,S0022IH2,S00622H1,S00BS4H1,S004BTH2,C000S5H2
## 3 S007DDH2,S007G7H4,C0066PH1,S00C2FH1,S00390H1,S001MJH1,S001S7H2,S0022IH2,S00622H1,S00BS4H1,S004BTH2
## 4 C0066PH1,S00C2FH1
## 5 C0066PH1,S00C2FH1
## 6 C000S5H1Con las coordenadas anteriores se espera conocer si existe mayor enriquecimiento de promotores o enhancers asociados con los resultados del experimento de CHIP-Seq de H3K4me1.
CHiP-seq histona H3k4me1
El proyecto Blueprint en epigenómica fue una iniciativa internacional que tuvo como objetivo mapear y caracterizar los perfiles epigenómicos de diferentes tipos celulares en la sangre.
El proyecto Blueprint utilizó diferentes tecnologías de secuenciación y análisis para estudiar la regulación génica a nivel epigenético. Estos incluyen la secuenciación de bisulfato de sodio, la cromatina inmunoprecipitación seguida de secuenciación (ChIP-seq), la detección de metilación del ADN y la secuenciación de ARN.
Para nuestro análisis se tomó un experimento de CHiP-seq de un neutrófilo. Los datos se encuentran en http://dcc.blueprint-epigenome.eu/#/experiments/ERX547970.
Procesamiento de los resultados del experimento CHIP-Seq
Se busca obtener los archivos BigWig y BED a partir de los datos del experimento.
Mapping
El mapeo se llevó a cabo utilizando bwa 0.7.7 para la
referencia GRCh38 del genoma humano.
bwa aln -q 5 grch38.fa input.fastq.gz > intermediate.sai ; bwa samse -r
"read group information" grch38.fa intermediate.sai input.fastq.gz |
samtools view -bS - > output.bam
# bwa aln: Comando bwa utilizado para alinear los datos de lectura a
#una referencia de genoma.
# -q 5: Se establece un umbral de calidad de base de 5 para las
#puntuaciones de calidad de llamada de base(bases con calidad menor
#a 5 no se alinearán).
# grch38.fa: Archivo de referencia.
# input.fastq.gz: Archivo de datos de lectura de entrada.
# > intermediate.sai: Este operador ">" redirige la salida del comando a
#un archivo llamado "intermediate.sai". El cual es una salida intermedia
#que se utilizará más adelante en la alineación.
# ";": Separador de comandos.
# bwa samse: comando BWA utilizado para generar un archivo SAM a
#partir del archivo SAI (lecturas single-end").
# -r "read group information": Agrega "información del grupo de lectura"
#al archivo SAM generado (útil en análisis posteriores).
#| samtools view -bS - > output.bam: Este comando redirige la salida
#del comando BWA samse a otro comando, samtools view.
# samtools view: convierte el archivo SAM generado por bwa samse a un
#archivo BAM.
# -bS: indica que la entrada es un archivo SAM y que la salida debe ser un
#archivo BAM.Los archivos BAM se clasifican y marcan como duplicados usando picard:
java -Xmx2048m -jar picard/SortSam.jar INPUT=input.bam
OUTPUT=output.bam SORT_ORDER=coordinate VALIDATION_STRINGENCY=SILENT
# java: comando para ejecutar la herramienta Picard, que se escribe
#en lenguaje de programación Java.
# -Xmx2048m: Establece el tamaño máximo de la memoria disponible para
#la herramienta Picard (2048 megabytes o 2 gigabytes).
# -jar picard/SortSam.jar: Especifica el archivo JAR (archivo de Java) que
#se utilizará para la herramienta Picard, en este caso SortSam.jar
#que se encuentra en el directorio "picard".
# INPUT=input.bam: Especifica el archivo de entrada que se utilizará
#para ordenar.
# OUTPUT=output.bam: Especifica el nombre del archivo de salida que se
#generará.
# SORT_ORDER=coordinate: Especifica el criterio de ordenación utilizado
#para ordenar las lecturas en el archivo de salida, en este caso por
#coordenadas de mapeo.
# VALIDATION_STRINGENCY=SILENT: Especifica el nivel de restricción que se
#utilizará para validar los registros de lectura en el archivo de
#entrada. En este caso "SILENT" permitirá la omisión de algunos registros
#que no cumplen con los requisitos de formato del archivo SAM/BAM.
java -Xmx2048m -jar picard/MarkDuplicates.jar INPUT=input.bam OUTPUT=output.bam METRICS_FILE=output.dup_metrics REMOVE_DUPLICATES=false ASSUME_SORTED=true VALIDATION_STRINGENCY=SILENT
# En este caso el archivo de java "MarkDuplicates.jar" se encuentra en el
#directorio picard.
# METRICS_FILE=output.dup_metrics: Indica el archivo de salida que
#contendrá las métricas generadas por el programa. En este caso
#"output.dup_metrics".
# REMOVE_DUPLICATES=false: Indica si se deben eliminar los duplicados o no.
#En este caso no.
# ASSUME_SORTED=true: Indica si se debe asumir que el archivo de entrada ya
#está ordenado. En este caso si.Filtering
A continuación, se filtró el archivo bam de salida para
eliminar las lecturas no asignadas y las lecturas con una calidad de
asignación inferior a 5:
samtools view -b -F 4 -q 5 input.bam > intermediate.output.bam
# -b: Especifica que la salida del comando será un archivo BAM
# -F 4: Deben excluirse los registros que tienen la flag "4" en el campo
#FLAG de la línea de registro SAM, o sea los que no se han alineado
#correctamente o no tienen alineaciones.A continuación, se filtró el archivo bam de salida
intermedia para eliminar PCR o lecturas duplicadas ópticas:
samtools view -b -F 1024 intermediate.output.bam > output.bam
# -F 1024: Excluye en este caso la flag "1024" de la línea de registro
#SAM. Se refiere a los registros marcados como duplicados por Picard.Modeling Fragment Size
El tamaño del fragmento se modela utilizando el script PhantomPeakQualTools R:
run_spp.R -c=output.bam -rf -out=params.out
# "run_spp.R": Es un script utilizado para ejecutar la herramienta de
#análisis de picos de phantom (SPP) en
#archivos de formato BAM.
# "-c": Especifica la ruta al archivo "BAM" (previamente obtenido) de
#entrada que se analizará.
# "-rf": Indica que se debe realizar una evaluación de la calidad de la
#lectura fantasma, necesario para determinar la escala y el número de
#picos a estimar.
# "-out": Especifica el nombre del archivo de salida que contendrá los
#resultados del análisis.Un script en R equivalente a run_spp.R
podría verse de la siguiente manera:
library(SPP)
library(PhantomPeakQualTools)
# Cargar "PhantomPeakQualTools" y "SPP".
input_bam <- "output.bam"
# Lee el archivo "BAM" filtrado anteriormente.
params_out <- "params.out"
# Archivo de salida que contiene los parámetros estimados
#del tamaño de fragmento.
spp_peak <- run_spp(input_bam, outputDir = ".", inputFile = input_bam,
treat = NULL, control = NULL, bin = 200,
blacklist = NULL, chrList = NULL, enrichment_thresh = 3,
norm_method = "RPM", filterDuplicates = TRUE, shift_len = 150,
cex_thresh = 0.9, corr_read_len = FALSE, correctForBigLen = TRUE,
min_len = 50, max_len = 1000, verbose = TRUE)
# "run_spp": Función en R que estima el número de picos de un
#experimento de secuenciación de ChIP-seq o CLIP-seq.
# "outputDir": Especifica el directorio de salida donde se guardarán los
#resultados. En este caso, el actual (".").
# "inputFile=input_bam": El archivo de entrada a utilizar.
# "treat" y "control" ="NULL" : se asume que no hay control y que todos
#los datos son tratables.
# "bin": Tamaño de la ventana que se utilizará para mapear los fragmentos.
# "blacklist": Archivo que contiene una lista de regiones que deben
#excluirse del análisis.
# "chrList": Archivo que contiene una lista de cromosomas que se
#utilizarán en el análisis.
# "enrichment_thresh": Umbral de enriquecimiento a utilizar para
#determinar el tamaño del fragmento.
# "norm_method": Método de normalización utilizado para los datos,
#en este caso, "RPM".
# "filterDuplicates": Indica si se deben filtrar los duplicados.
# "shift_len": Número de bases para el desplazamiento de fragmentos.
# ""cex_thresh": Umbral de coeficiente de correlación a utilizar para
#determinar la calidad del desplazamiento del fragmento.
# "corr_read_len": Indica si se debe corregir la longitud del
#fragmento en la correlación.
# "correctForBigLen": Indica si se debe corregir la longitud del fragmento
#para fragmentos grandes.
# "min_len" y "max_len": Tamaño mínimo y máximo de fragmento a incluir en
#el análisis.
# "verbose": Indica si se deben imprimir mensajes de progreso detallados.
write.table(as.data.frame(t(spp_peak)), file = params_out, quote = FALSE, sep = "\t", row.names = FALSE)
# "write.table": Escribe objetos de datos en un archivo de texto. En este
#caso "fragments$FragmentLengths" que contiene los tamaños de
#fragmentos obtenidos, en el archivo "fragment_sizes.txt".
# "as.data.frame": Convierte el objeto de salida de la función "run_spp"
#en un marco de datos (data frame) de R.
# "t(spp_peak)": Realiza la transposición(intercambia sus filas por sus
#columnas) del objeto "spp_peak" .
# "file": Indica el nombre del archivo de salida.
# "sep": Especifica el separador de columna utilizado en el archivo de
#salida (\t=tabulador).
# "quote = FALSE": Evita que se utilicen comillas alrededor de los
#valores en el archivo de salida.
# "row.names": Evita que se incluyan los nombres de las filas en el
#archivo de salida.Peak Calling
MACS2 se usa para llamadas máximas con el tamaño de
fragmento predicho por PhantomPeakQualTools. Utilizan tanto
el método de carrera estándar como la -broad flag
dependiendo de la marca en cuestión.
macs2 callpeak -t chip.bam -n a_sensible_name --gsize hs -c input.bam --nomodel --shiftsize=half_fragment_size --broad
# "-t chip.bam": Es el archivo BAM de las lecturas de ChIP-seq.
# "-n a_sensible_name": Nombre significativo para la salida de los
#picos llamados.
# "--gsize hs": Indica el tamaño del genoma de la especie en estudio.
#En este caso del genoma humano.
# "-c input.bam": Archivo "BAM" de las lecturas de control.
# "--nomodel": Indica que no se ajustará un modelo de tamaño de fragmento
#específico para los datos de ChIP-seq, sino que se utilizará una
#aproximación constante para el tamaño del fragmento.
# "--shiftsize=half_fragment_size": Tamaño del corrimiento para el
#análisis de los picos, en este caso la mitad del tamaño del fragmento
#obtenido.
# "--broad": especifica que se buscan picos anchos en lugar de picos
#estrechos.
macs2 callpeak -t chip.bam -n a_sensible_name --gsize hs -c input.bam --nomodel --shiftsize=half_fragment_size
# En este caso se buscan picos estrechos.Las marcas en las que se utilizó -broad son:
H3K27me3
H3K36me3
H3K9me3
H3K4me1Las marcas donde se omite -broad son:
H3K27ac
H3K4me3
H3K9/14ac
H2A.ZacBedGraph
Los gráficos de señal donde el signo se traza en el eje vertical y la
posición a lo largo del genoma se traza en el eje horizontal, se
producen usando align2RawSignal y el tamaño del fragmento
predicho por PhantomPeakQualTools. Se utilizan archivos
fasta y umap específicos de sexo.
align2rawsignal -i=chip.bam -of=bg -o=chip.bg -l fragment_size -s=/path/to/fasta_files -u=/path/to/umap_files
# "align2rawsignal": Genera un archivo de señal(usa para "wiggle plots")
#con información sobre la cantidad de secuencias mapeadas en cada
#posición del genoma a partir de un archivo "BAM".
# "-i=chip.bam": Archivo BAM de entrada que contiene las
#lecturas de ChIP-seq.
# "-of=bg": Especifica el formato de salida de la señal cruda. En este
#caso "bg" o sea un archivo de señal "bedGraph".
# "-o=chip.bg": Nombre del archivo de salida que se generará.
# "-l fragment_size": Tamaño de fragmento utilizado en el análisis de
#ChIP-seq, obtienido mediante "PhantomPeakQualTools".
# "-s=/path/to/fasta_files": La ruta al archivo fasta correspondiente al
#genoma de referencia utilizado en el ChIP-seq.
# "-u=/path/to/umap_files": Especifica la ruta al archivo UMAP
#correspondiente al genoma de referencia utilizado en el ChIP-seq.Del código anterior se obtiene un archivo “bedGraph”(bg) que se convierte en los tipos de archivos que utilizamos; BED y BigWig (bw) de la siguiente manera:
bedGraph a BigWig:
1- Descarga el programa
bedGraphToBigWig del directorio de utilidades binarias (http://hgdownload.soe.ucsc.edu/admin/exe/).
2- Usa el script fetchChromSizes del
mismo directorio para crear el archivo chrom.sizes para la base
de datos UCSC con la que está trabajando. Si el ensamblado
genNom está alojado en UCSC, chrom.sizes puede ser
una URL;
http://hgdownload.soe.ucsc.edu/goldenPath/genNom/bigZips/genNom.chrom.sizes.
3- Use bedGraphToBigWig de la siguiente
manera para crear un archivo bigWig a partir de un bedGraph(sin
comprimir):
bedGraphToBigWig in.bedGraph chrom.sizes myBigWig.bw4- Mueve el archivo bigWig a una ubicación http, https o ftp accesible desde la web.
5- Pegue la URL en el formulario de entrada de pista personalizada o construya una pista personalizada(más información->https://genome.ucsc.edu/goldenPath/help/hgTracksHelp.html#CustomTracks) utilizando una sola línea de pista(más información->https://genome.ucsc.edu/goldenPath/help/hgTracksHelp.html#TRACK).
6- Pegue la línea de seguimiento personalizado en el cuadro de texto de la página de administración de seguimiento personalizado(https://genome.ucsc.edu/cgi-bin/hgCustom).
BigWig a Bed:
1- Descarga el archivo bigBedToBed desde el directorio de utilidades binarias de UCSC Genome Browser (http://hgdownload.soe.ucsc.edu/admin/exe/).
2- Ejecuta el comando siguiente en la terminal:
bigBedToBed in.bigWig out.bed
# "in.bigWig": Nombre del archivo bigWig de entrada.
# "out.bed": Nombre del archivo BED de salida. Se creará en el
#directorio actual.Los resultados de la llamada de picos están dados en la siguiente tabla con formato .bed:
H3K4me1_neutrophil_2 <- read.table("H3K4me1_neutrophil_2.bed.gz", header = F, sep = "\t")
head(H3K4me1_neutrophil_2)## V1 V2 V3 V4 V5 V6 V7 V8
## 1 chr1 190784 192075 3648.macs2_peak_call_peak_1 118 . 7.78170 13.70875
## 2 chr1 264665 264911 3648.macs2_peak_call_peak_2 31 . 4.18846 4.74141
## 3 chr1 267976 268542 3648.macs2_peak_call_peak_3 18 . 3.43295 3.37943
## 4 chr1 270477 271574 3648.macs2_peak_call_peak_4 52 . 5.11515 6.92670
## 5 chr1 272744 273847 3648.macs2_peak_call_peak_5 18 . 3.40187 3.28042
## 6 chr1 778023 779737 3648.macs2_peak_call_peak_6 29 . 4.00693 4.50536
## V9
## 1 11.89800
## 2 3.18868
## 3 1.89258
## 4 5.29479
## 5 1.80105
## 6 2.96993Información de los Bed files:
1- chrom
2- start
3- end
4- name
5- score (fold enrichment * 10, rounded down to integer
value)
6- strand
7- -log10(qvalue)
8- -log10(pvalue)
9- fold_enrichment
Enriquecimiento de promotores vs enhancers
Regiones promotoras y enhancers
Como ya se ha descrito, en el archivo “ActivePromoterEnhancerLinks.tsv”, encontramos los promotores y los enhancers activos en las celulas sanguineas, por lo que se extraen las coordenadas de ambos grupos de la siguiente forma:
library(GenomicRanges)
# obtener las coordenadas de los promotores, fragmentos
GR_promotores <- makeGRangesFromDataFrame(interacciones, keep.extra.columns = T,
start.field = "baitSt", end.field = "baitEnd",
seqnames.field = "baitChr", strand.field = "*")
GR_fragmentos <- makeGRangesFromDataFrame(interacciones, keep.extra.columns = T,
start.field = "oeSt", end.field = "oeEnd",
seqnames.field = "oeChr", strand.field = "*")Se exporta a formato bed para hacer una conversion en UCSC a la referencia hg38, ya que el experimento fue alineado con la referencia hg19.
export.bed(GR_promotores,con='promotores.bed')
export.bed(GR_enhancers,con='enhancers.bed')En la pagina https://genome.ucsc.edu/cgi-bin/hgLiftOver se hace la conversion.
Se eliminan las regiones repetidas en la terminal para obtener las regiones únicas con el siguiente comando.
awk '!x[$0]++' promotores_hg38.bed > promotores_unicos.bed
awk '!x[$0]++' enhancers_hg38.bed > enhancers_unicos.bedVisualización de datos
En el archivo “H3K4me1_neutrophil_2.bed.gz” se encuentran los picos de la señal del experimento chip-seq de la histona, usamos este porque ya esta procesado. Se pueden visualizar todos estos archivos en IGV:
Intersecciones
Con estos archivos se buscan las intersecciones de los promotores y los enhancers con la señal de los picos de histona:
# regiones de interseccion:
# Se obtienen las coordenadas de promotores y enhancers completas caundo tienen la modificación
bedtools intersect -wa -a promotores_unicos.bed -b H3K4me1_neutrophil_2.bed.gz | awk '!x[$0]++' > promotores_activos_c.bed
bedtools intersect -wa -a enhancers_unicos.bed -b H3K4me1_neutrophil_2.bed.gz | awk '!x[$0]++' > enhancers_activos_c.bedEs necesario usar la opcion -wa para tomar la region de
promotores y enhancers completa y así poder ver la distribución de la
modificación de histona a lo largo de las secuencias.
computeMatrix scale-regions -S H3K4me1_neutrophil_2.bw -R promotores_activos_c.bed --beforeRegionStartLength 4000 --regionBodyLength 4000 --afterRegionStartLength 4000 --skipZeros -o matrix_promotores_H3K4me1.tab.gz
plotHeatmap -m matrix_promotores_H3K4me1.tab.gz -out H3K4me1_promotores.png --regionsLabel enhancers
computeMatrix scale-regions -S H3K4me1_neutrophil_2.bw -R enhancers_activos_c.bed --beforeRegionStartLength 4000 --regionBodyLength 4000 --afterRegionStartLength 4000 --skipZeros -o matrix_enhancers_H3K4me1.tab.gz
plotHeatmap -m matrix_enhancers_H3K4me1.tab.gz -out H3K4me1_enhancers.png --regionsLabel enhancers --plotTitle H3K4me1.enhancers
En el heatmap de promotores observamos un enriquecimiento a lo largo de la secuencia, sin embargo podemos notar que este enriquecimiento se encuentra más remarcado en los enhancers pues a los alrededores de dicha región dejan de estar las modificaciones.
Promotores/Enhancers vs H3k4me1
Con el número de enhancers y promotores que están activos se puede crear una tabla de contingencia para hacer un chi-test y, en primer lugar probar si existe una relación entre promotores/enhancers y las marcas de histona.
a <- 7656 # número de promotores asociados con la histona
b <- 8017-7656 # número de promotores no asociados con la histona
c <- 10757 # número de enhancers asociados con la histona
d <- 25225-10757 # número de enhancers no asociados con la histona
# Construir la tabla de contingencia como una matriz
tabla_contingencia <- matrix(c(a, b, c, d), nrow = 2, byrow = TRUE,
dimnames = list(c("Promotores", "Enhancers"),
c("Histona", "No Histona")))
tabla_contingencia## Histona No Histona
## Promotores 7656 361
## Enhancers 10757 14468La prueba chi-cuadrado es útil en este caso porque se utiliza para evaluar la relación entre dos variables categóricas, que es lo que se está analizando en la tabla de contingencia de promotores y enhancers asociados y no asociados con la histona.
En este contexto, se está evaluando si hay una diferencia significativa en la proporción de promotores y enhancers que están asociados con la histona en comparación con los que no lo están. La hipótesis nula de la prueba chi-cuadrado es que no hay relación entre las dos variables categóricas (en este caso, la histona y los promotores/enhancers). Si se rechaza la hipótesis nula, se concluye que hay evidencia suficiente para sugerir que hay una relación significativa entre las variables.
Para hacer la prueba se usa la función chisq.test en
R.
chisq.test(tabla_contingencia)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tabla_contingencia
## X-squared = 6875.4, df = 1, p-value < 2.2e-16Dado que el valor del p-valor es muy pequeño, se puede concluir que la relación entre los promotores/enhancers y la histona es muy significativa y no es probable que sea el resultado de una casualidad.
Promotores vs Enhancers
Sabiendo que existe relación entre ambas variables, para comparar el número de promotores y enhancers, una prueba adecuada es la prueba de chi-cuadrado de homogeneidad. Esta prueba se utiliza para determinar si hay una diferencia significativa en la proporción de individuos en diferentes grupos.
Para realizar la prueba, se tiene que contar el número de promotores y enhancers en los datos y compararlos. Luego, usar la prueba de chi-cuadrado de homogeneidad para determinar si hay una diferencia significativa entre los grupos en términos de su proporción.
total_promotores <- 7656
total_enhancers <- 10757
tabla_contingencia <- matrix(c(total_promotores, total_enhancers), nrow = 2, byrow = TRUE)
chisq.test(tabla_contingencia)##
## Chi-squared test for given probabilities
##
## data: tabla_contingencia
## X-squared = 522.25, df = 1, p-value < 2.2e-16Dado que el p-value es muy pequeño, se puede concluir que hay una diferencia significativa entre los grupos en términos de su proporción y como hay mas enhancers podemos concluir que estos son los mas enriquecidos.
Enriquecimiento de factores de transcripción
Sabemos que existen promotores y enhancers que se encuentran activos,
ahora necesitamos saber que factores de transcripción (TF) que se unen a
estos sitios. Para ello se utiliza la función matrix-scan
de rsat.
El objetivo de esto es analizar los TFs que se unen a lo largo de todos los promotores y enhancers encontrados; sin embargo esto es una tarea que necesita mucho tiempo de ejecución. Es por esto que hemos decidido únicamente buscar los TFs que se encuentran en el cromosoma 22.
Nota: Para el cómputo de enriquecimiento en enhancers se dividió el archivo en seis partes para disminuir aún más el tiempo de cómputo.
Para realizar el análisis de enriquecimiento se necesitan varios archivos.
Matrices de Posicion Frecuencia (PFMs) de la colección de JASPAR.
La base de datos JASPAR contiene matrices de PFM para una amplia variedad de factores de transcripción en diferentes especies, incluyendo humanos. Las matrices JASPAR se crean a partir de datos experimentales de sitios de unión de factores de transcripción y se utilizan para predecir la ubicación y la fuerza de la unión de un factor de transcripción a un sitio de ADN específico.
Existen dos colecciones principales: la colección no redundante y la redundante. La colección no redundante de JASPAR se refiere a un conjunto de factores de transcripción que han sido agrupados de tal manera que no hay superposición entre ellos en términos de secuencias de unión al ADN.
Por otro lado la colección redundante de JASPAR contiene múltiples secuencias de unión al ADN para cada factor de transcripción, lo que significa que puede haber cierta superposición entre ellas.
Ambas colecciones se pueden obtener en estos links:
Secuencias de promotores y enhancers
Cada matriz debe escanearse sobre una secuencia para determinar la probabilidad de que esa secuencia sea un sitio de unión a TF. Para obtener estas secuencias se usó el siguiente código.
# obtener solo el cromosoma 22
awk '$1 == "chr22"' promotores_activos.bed > promotores_activos_chr1.bed
awk '$1 == "chr22"' enhancers_activos.bed > enhancers_activos_chr1.bed
# con la erramienta bedtools obtener archivos fasta
bedtools getfasta -fi hg38.fa -bed promotores_activos_chr22.bed -fo promotores_activos_chr22.fa
bedtools getfasta -fi hg38.fa -bed enhancers_activos_chr22.bed -fo enhancers_activos_chr22.faTambién se usará una secuencia permutada como control la cual se generó de esta manera:
random-seq -l 1000 -n 10 -format fasta -expfreq 2nt_upstream-noorf_Homo_sapiens_GRCh37-ovlp-1str.freq.gz -o secuenciaspermutadas.fa
# -l : tamaño de la secuencia
# -n : número de secuencias
# -expfrec modelo de frecuencia de fondoLa secuencia permutada tiene el mismo modelo de fondo para tener la misma proporción de nucleotidos que las secuencias originales.
Sin embargo, lo ideal es permutar las secuencias originales con este código:
random-seq -template_format len -i promotores_activos_chr22.fa -format fasta -expfreq 2nt_upstream-noorf_Homo_sapiens_GRCh37-ovlp-1str.freq.gz -o secuenciaspermutadas.faFrecuencia de fondo
La frecuencia de fondo es la frecuencia esperada de ocurrencia de cada base en una secuencia de ADN. Esta frecuencia se basa en la distribución de bases en el genoma de referencia.
La frecuencia de fondo se utiliza en muchos análisis de secuencias de ADN para estimar la probabilidad de que una secuencia determinada se encuentre en una región no codificante o no reguladora del genoma. Al comparar la frecuencia observada de ocurrencia de una secuencia con la frecuencia de fondo esperada, se puede determinar si la secuencia está enriquecida o empobrecida en una muestra determinada.
En el código que se usa más adelante es posible usar un modelo de fondo dado un contexto génico, en nuestro caso usamos este:
En él se observa la probabilidad de obtener un cierto par de nucleótidos en el genoma hg38.
Es posible usar otros modelos tomando en cuenta más nucleótidos e incluso usar cadenas de Markov.
Si no queremos usar un archivo de frecuencia de fondo, este puede calcularse automáticamente dado el archivo fasta.
Umbral del p-valor
Al realizar el escaneo es muy probable obtener falsos positivos, es
decir, predecir un sitio de unión a TFs cuando no lo es. Esto debido a
que existen matrices que son más pequeñas y por tanto tienen más
probabilidad de aparecer de forma aleatorea en el genoma o en el caso de
otras matrices pueden ser muy variables. Es por eso que se necesita un
p-valor mínimo para disminuir este error. Para ello usamos la función
matrix-distrib de rsat:
matrix-distrib -v 1 -top 100 -m JASPAR2022_CORE_vertebrates_non-redundant_pfms_jaspar.txt -pseudo 1 -decimals 1 -matrix_format jaspar -bgfile 2nt_upstream-noorf_Homo_sapiens_GRCh37-ovlp-1str.freq.gz -bg_pseudo 0.01 -o md_data.tab
XYgraph -v 1 -title1 "Distribution of weights (log scale)" -title2 "Score probability and P-value" -pointsize 2 -lines -symbols -legend -bg white -xcol 1 -xleg1 "weight" -xmin -10 -xmax 20 -xsize 500 -ycol 2,4 -yleg1 "frequency (log 10 scale)" -ysize 500 -ylog 10 -i md_data.tab -format png -o XYgraph.png -htmap 
En esta gráfica podemos ver en el eje X la distribución de los pesos de las primeras 100 matrices de la colección no redundante, esta nos indica si la secuencia escaneada corresponde a un factor de transcripción que corresponde a la matriz o es una secuencia dada solamente por el contexto genómico. Valores mayores a 5 son buenos para decidir que sí es un sitio de unión.
En verde podemos ver la distribución del p-valor en escala logarítmica, esta nos ayuda a definir el p-valor necesario para disminuir la cantidad de falsos positivos.
Un p-valor de 1e-4 o 1e-5 podría ser suficiente para disminuir el error de las matrices, esto significa que cada 10000 bases encontraríamos una secuencia a la que no se une un factor de transcripción pero diríamos que sí.
Matrix-scan
Con lo anterior podemos escanear las matrices usando el siguiente código:
matrix-scan -v 1 -matrix_format jaspar -m archivo_de_matrices.txt -pseudo 1 -decimals 1 -2str -origin end -bgfile archivo_de_frecuencia.freq -return distrib -return occ_proba -sort_distrib -uth occ_pval 1e-4 -lth occ_sig 0 -uth occ_sig_rank 1 -lth score 5 -i archivo_fasta.fa -seq_format fasta -n skip -o lista_de_ptomotores.tab
# -v : activa el modo detallado de mensajes de registro para el programa.
# -matrix_format : especifica el formato de la matriz
# -m : indica el archivo de texto que contiene las matrices de puntuación utilizadas para el análisis.
# -pseudo 1 : agrega una pseudocuenta a las frecuencias de la matriz para evitar problemas de cálculo cuando las frecuencias son muy pequeñas.
# -decimals 1: especifica el número de decimales a utilizar en la salida.
# -2str : convierte las secuencias de ADN en ambas hebras, esto por que hay factores de transcripción que se unen en ambas cadenas.
# -origin end : especifica que los sitios se identificarán a partir del final de la secuencia, esto hace que las posiciones de coincidencia sean valores negativos, proporcionando la distancia entre la coincidencia y el final de la secuencia.
# -bgfile: especifica un archivo que contiene las frecuencias de fondo para cada base del ADN.
#-return distrib : devuelve la distribución de puntuación
# -return occ_proba": devuelve la probabilidad de ocurrencia de cada sitio. Esto para el enriquecimiento de TF en todo el conjunto de secuencias de entrada
# -sort_distrib : ordena la distribución de puntuación.
# -uth occ_pval : establece un umbral de p-valor
# -lth occ_sig : establece un umbral de valor para la significancia de la ocurrencia de los sitios
# -uth occ_sig_rank : establece un umbral para la clasificación de significancia de los sitios, es decir, devuelve el mejor en cada secuencia
# -lth score 5 : establece un umbral para la puntuación de los sitios.
# -i : especifica el archivo que contiene las secuencias de ADN a analizar.
# -seq_format : especifica el formato de las secuencias de ADN.
# -n skip : Las regiones que contienen N se omiten.
# -o : especifica el nombre del archivo de salida que contendrá la lista de sitios de unión de factores de transcripción identificados y ordenados con mayor probabilidad y enriquecimiento¿Por qué me conviene usar la colección no redundante?
Tabla de enriquecimiento de TF en promotores con la colección no redundante
Tabla de enriquecimiento de TF en promotores con la colección redundante
Al comparar ambas tablas no es posible argumentar a simple vista la colección más conveniente para utilizar, sin embargo nos conviene usar una matriz no redundante de JASPAR porque contiene sólo una matriz por cada familia de factores de transcripción, lo cual reduce la complejidad de los datos, lo que facilita la interpretación y la identificación de patrones de unión específicos.
Al eliminar matrices redundantes, se evita la sobreestimación de la frecuencia de los sitios de unión de los factores de transcripción en una secuencia dada, que puede llevar a la identificación de sitios de unión falsos.
Promotores vs Enhancers con colección no redundante
Los TFs que comparten promotores y enhancers con matrices no redundates, se muestran a continuación:
Esta tabla muestra el enriquecimiento de los factores de transcripción y sus scores, lo cual indica qué tan asociados están dichos TFs con los enhancers.
Esta tabla muestra los factores de transcripción asociados a los promotores.
Como se observa, ambas tablas comparten TFs en común, tales
como:
1) Arnt
2) Ahr::Arnt
3) Ddit3::Cebpa
4) MAX::MYC
5) RORA
6) RREB1
7) ZNF354C
8) CTCF
9) EWSR1-FLI1
10) PLAG1
11) Zfx
12) FOS
13) FOSL2
14) Gfi1B
Entre muchos más TFs, que podrían indicar la interacción entre los promotores y los enhacers.
¿Por qué usar secuencias permutadas me sirve como control?
El uso de secuencias permutadas como control ayudaría a evaluar qué tan significativamente enriquecida está la secuencia en sitios de unión de factores de trancripción, aleatorizando la posición de las bases de la secuencia original.
Al comparar la frecuencia de unión del TF en la secuencia original con la secuencia permutada, se puede determinar si la unión es específica de la secuencia original o si es sólo un resultado de una distribución aleatoria de sitios de unión en la secuencia.
Si la frecuencia de unión del factor de transcripción en la secuencia original es significativamente mayor que la frecuencia de unión en las secuencias permutadas, entonces se puede concluir que la secuencia original contiene un sitio de unión específico para ese factor de transcripción.
Podemos observar lo dicho anteriormente en el siguiente ejemplo:
Esta imagen corresponde a una tabla de enriquecimiento de los factores
de transcripción y sus scores que muestran qué tan asociados están con
las secuencias originales (enhancers) con matrices de JASPAR no
redundantes.
Por otra lado, esta imagen contiene el control que son secuencias permutadas con matrices de JASPAR no redundantes.
Como se puede ver, ambas tablas comparten algunos factores de transcripción, tales como Arnt y RORA. Esto quiere decir que dichos TFs son únicamente resultado de la distribución aleatoria de las secuencias de las matrices de JASPAR y no un sitio de unión específico.
Conclusión
Hemos detectado un mayor enriquecimiento de enhancers asociados con la modificación H3K4me1 en comparación con los promotores, lo que concuerda con los datos esperados respecto al artículo de Javierre et al., 2016.
También encontramos los sitios de unión a factores de transcripción en el cromosoma 22. Sin embargo, hace falta un análisis más profundo para determinar si estos sitios de unión reflejan algún tipo de función en relación con su tipo celular.
Por último, destacamos que el uso de la colección no redundante puede ayudar a simplificar la interpretación de los resultados y que el uso de secuencias permutadas reduce el número de falsos positivos.
Bibliografía
- Biola M. Javierre, Oliver S. Burren, Steven P. Wilder, …, Chris Wallace, Mikhail Spivakov, Peter Fraser (2016). 2) Lineage-Specific Genome Architecture Links Enhancers and Non-coding Disease Variants to Target Gene Promoters. Cell 167, 1369-1384.
- OpenAI. (2021). ChatGPT. https://openai.com/blog/chat-gpt-3-language-models/
- RSAT Metazoa: http://rsat.sb-roscoff.fr/tutorials.php
- Rsat matrix-scan output file format documentation: http://www.ra.cs.uni-tuebingen.de/software/ModuleMaster/doc/html/node6.html
- bedtools 2.4.20: https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html
- deeptools 2.5.3: https://deeptools.readthedocs.io/en/develop/
- JASPAR 2022: https://jaspar.genereg.net/
- matrix-scan: http://rsat.sb-roscoff.fr/help.matrix-scan.html
- UCSC Genome Bioinformatics. (s.f.). Table Browser. Recuperado em de abril de 2023, de http://hgdownload.soe.ucsc.edu/admin/exe/
- Blueprint Data Analysis Center. (s.f.). Blueprint Epigenome Project. Recuperado en abril de 2023, de http://dcc.blueprint-epigenome.eu/#/experiments/ERX547970
- R 4.2.2: https://cran.r-project.org/